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Abstract 

This work is based upon a coupled, lattice-based continuum formu- 
lation that was previously applied to problems involving strong coupling 
between mechanics and mass transport; e.g. diffusional creep and elec- 
tromigration Here we discuss an enhancement of this formulation 

to account for migrating grain boundaries. The level set method is used 
to model grain-boundary migration in an Eulerian framework where a 
grain boundary is represented as the zero level set of an evolving higher- 
dimensional function. This approach can easily be generalized to model 
other problems involving migrating interfaces; e.g. void evolution and free- 
surface morphology evolution. The level-set equation is recast in a re- 
markably simple form which obviates the need for spatial stabilization 
techniques. This simplified level-set formulation makes use of velocity ex- 
tension and field re-initialization techniques. In addition, a least-squares 
smoothing technique is used to compute the local curvature of a grain 
boundary directly from the level-set field without resorting to higher- 
order interpolation. A notable feature is that the coupling between mass 
transport, mechanics and grain-boundary migration is fully accounted for. 
The complexities associated with this coupling are highlighted and the 
operator-split algorithm used to solve the coupled equations is described. 
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1 Introduction 



In a recent paper, Garikipati et al. pQ presented a coupled continuum field for- 
mulation for the interaction of electric effects, mechanical response and self- 
diffusion. Their approach drew upon earlier work by Larche and Cahn 
Nix and co-workers [SJ|H], Genin [7], Bower and Freund Xia et al. [5] — to 
name but a few. In Reference pQ, a review of these works and others was pre- 
sented, computational techniques were developed based on the finite element 
method, and several initial and boundary value problems involving diffusional 
creep (Nabarro-Herring and Coble creep) were solved. 

The formulation introduced in was extended to interdiffusion, with dopants 
in silicon as motivation, by Garikipati and Bassman In the current pa- 
per, we present a further extension of that same formulation to account for 
the interaction of grain-boundary migration with stress-driven self-diffusion and 
electromigration in polycrystalline solids. 

Different strategies used to model grain-boundary motion in a computational 
setting are described in the literature. Sun and Suo jlOl ITT] developed a two- 
dimensional finite element formulation, based on the idea that the energy dissi- 
pated during the motion of the boundary must equal the reduction in the free 
energy of the system. This formulation was used to study several problems in- 
cluding grain growth in a thin film and the competition between surface grooving 
and grain-boundary migration. A similar variational formulation was developed 
by Cocks and Gill ^1 Ej and used to model the evolution of a large network 
of grains in two dimensions. Another model was developed by Zhao et al. |14| . 
based on the variational formulation of Reitich and Soner |15] and using the 
level set method of Osher and Sethian |16j . 

The level set method is also used in the present paper as an interface-capturing 
technique. However, the goal of the current work is not merely to simulate grain- 
boundary motion, but to capture fully the interaction between this motion and 
other microscale phenomena that take place in pure polycrystalline materials, 
namely stress-mediated self-diffusion and electromigration. This distinguishes 
the current work from the existing literature. 

In the current work, we recast the level-set equation in a simpler form by assum- 
ing that the level-set function remains a signed distance to the migrating grain 
boundary (as with the original level-set equation, the use of an extensional veloc- 
ity field helps maintain this signed-distance function). Mourad et al. |2J tested 
the resulting level-set formulation extensively and concluded that it is both ac- 
curate and robust despite its remarkable simplicity (for some interface-evolution 
problems, this approach reduces the original level-set equation, a nonlinear hy- 
perbolic PDE, to an ODE that is almost trivial to solve). They conducted sev- 
eral numerical experiments to assess the ability of the simplified level-set scheme 
to capture the correct solution, particularly in the presence of discontinuities in 
the extensional velocity and/or in the gradient of the level-set function. They 
also examined the convergence properties of the method and its performance in a 
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variety of problems, including curvature flow and problems where the simplified 
level-set equation takes the form of a Hamilton- Jacobi equation with convex or 
non-convex Hamiltonian. Discretizations based on structured and unstructured 
finite-element meshes of bilinear quadrilateral and linear triangular elements 
were shown to perform equally well. They also found that sufficient accuracy 
is available through a standard Galerkin formulation without resorting to any 
stabilization or discontinuity-capturing [181 119| techniques. 

In addition, a variant of the simplified level-set formulation mentioned above 
was employed by Ji et al. [201 f° r representing the evolution of phase boundaries 
over unstructured finite-element meshes. Here, we treat a complex coupled 
problem of which grain-boundary migration is only one facet, in addition to 
mechanics, self-diffusion and electromigration. We apply this simplified level-set 
scheme to the problem of grain-boundary migration within this context, and we 
demonstrate its implementation as an integral part of the wider computational 
framework used to solve the coupled problem under consideration. 

The remainder of this paper is organized as follows. In Section [21 we summarize 
the formulation for coupled mass transport and mechanics. Then we examine 
the thermodynamics and kinetics of grain-boundary migration, and show how 
this phenomenon interacts with mass transport and mechanics in polycrystals. 
In Section[31 we formulate the grain-boundary migration problem using the level 
set method and we describe the computational methods used in the implemen- 
tation of this formulation. Numerical examples are presented in Section A 
summary is provided and conclusions are drawn in Section [3J 

2 The coupled formulation 

2.1 Thermodynamic basis 

The thermodynamics is posed in a continuum setting, with motivation provided 
by atomic processes. A schematic of the lattice is shown in Fig. ^ It shows 
atoms, vacancies, and a free surface; the latter is a source and sink for vacancies. 
A grain boundary could serve as a source or sink also. Free surfaces and grain 
boundaries are treated as regions of finite width, S s and 25 g b respectively. 

2.1.1 Internal energy density 

We consider crystalline materials in which the dangling bonds around a vacancy 
cause an inward relaxation of the surrounding lattice (see Fig.^l. The resulting 
vacancy relaxation strain can be expressed as 

e v = -±{l-fMC v -C%)l, (1) 
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Figure 1: Schematic rendering of a lattice with atoms and vacancies. 



where fi is the atomic volume, /fl (with < / < 1) is the volume of a vacancy, 
C„ is the vacancy concentration, C% q is the vacancy concentration at thermo- 
dynamic equilibrium under vanishing external stress and 1 is the second-order 
isotropic tensor. The creep strain resulting from the accumulation or deple- 
tion of atoms at a free surface or grain boundary with unit normal, n, can be 
expressed as 

e c = \e c {n® n). (2) 

A non-phenomcnological evolution equation for 9 C was derived and discussed in 
detail by Garikipati et al. £Q. The thermal strain is given by 

e th =a(T-T )l, (3) 

where T is the temperature, To is a reference temperature and a is the linear 
coefficient of thermal expansion. The elastic strain is obtained by subtracting 
these inelastic strain contributions from the total strain, e; i.e. 

e e =e-(e v +e c + e th ). (4) 

The stress is obtained from the elastic strain and the (generally anisotropic) 
fourth-order elasticity tensor, C, as 

er = C:£ e . (5) 

Given a state with local stress, er, the incremental elastic strain energy density is 
then given as u : Se e — (C : e e ) : Se e , where 5e e is the increment in elastic strain. 
With this background, the incremental internal energy density 5e, corresponding 
to a state {e, T, C v }, and increments Se and 5C V is 

5e = 5e ri (rj) + eJ v 5C v + (C : e e ) : 6e e , (6) 
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where r\ is the entropy density and e{ is the vacancy formation energy. The 
specific form of e n (r]), the entropic dependence of e, is not important to the 
development that follows. 



2.1.2 External work density 

The density of work done by external agents can be expressed as 

3 

Swext — cr : Se - qijj5C v + -(n • crn)fflSC v x- (7) 

Here, the first term of the right-hand side is the classical stress-power term. 
The second term accounts for the apparent work performed by the electrostatic 
potential, tp, during electromigration; in this phenomenological treatment, q is 
the apparent charge ascribed to each vacancy. The last term accounts for the 
work done against the stress when vacancies are created at a free surface or grain 
boundary. The numerical factor appearing in this term arises from geometrical 
considerations. Additionally, the requirement that this term be active only at 
sources and sinks is enforced using the indicator function, x, defined as 

, , J 1 if a; is in a surface- or grain-boundary region, , . 
WMJ = < o otherwise. [ > 



2.1.3 Entropy density 

The total entropy density is given by 

V = Vvib(T) - k 



(9) 



Since the formulation is isothermal, the specific form of the vibrational term, 
Vvib(T), is unimportant. The second term is the entropy density due to mixing 
(see Kittel and Kroemer [5T] for details), k is the Boltzmann constant, C a is the 
concentration of atoms, and C s is the lattice site concentration. Assuming that 
C s remains fixed in any material volume, i.e. SC S = SC a + 5C V = 0, the incre- 
mental entropy density corresponding to an increment in vacancy concentration, 
SC V , at given stress and temperature can be expressed as 

Sr, = -k\og(^)sC v . (10) 



2.1.4 The Gibbs free energy density 

At the given state of stress and temperature, the incremental Gibbs free energy 
density corresponding to increments Se and SC V is defined as 

Sg = Se - Suiext ~ TSr\. (11) 
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With this, the constitutive relations can be obtained in a systematic fashion as 
outlined in the following section. 

Remark 1 Since free surfaces and grain boundaries are considered, there are 
accompanying surface and grain boundary energies, 7 S and y go . These energies 
are taken to be independent of the strain and vacancy concentration for this 
formulation, and therefore do not appear in the incremental Gibbs free energy 
density. 

2.2 Constitutive relations 

Anisotropic elasticity is assumed and the relation between the stress, <x, and 
the elastic strain, e e , is given by JSJ). The relation between the current density, 
z, and the electric potential, ip, is given by Ohm's law: 

where p is the electric resistivity. 

2.2.1 The chemical potential of vacancies 

The vacancy chemical potential is defined in the usual fashion |22|: 

fi v 5C v = n e v q 6C v + 6g, (13) 
where [i eq is a constant reference potential. Applying l|13fl to (|j) Hllf> gives 

fi v = n e v q + el + (C : e e ) : ±(1- f)Sll- ^(n- an) fn X +qi>+ kT log f ^ . (14) 

The coupling with mechanics is evident through the strain- and stress-dependent 
terms. 

2.2.2 The vacancy flux 

The vacancy flux is obtained from (|14fl via the relation 

where D v is the vacancy diffusivity. This relation has been derived from an 
atomic basis by Bardeen [231 - The term D v C v /kT is the mobility of a vacancy 
and — V/x„, the force acting on it, is the driving force for mass transport. 
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2.3 Governing equations 



The constitutive relations established above are incorporated in balance laws 
for mechanics, electric flow and mass transport, leading to a coupled system of 
governing differential equations. 

2.3.1 Mechanics 

Neglecting dynamic effects and body forces, the mechanics problem is governed 
by the quasistatic equilibrium equation and appropriate boundary conditions: 

V<x = 0, (16a) 
u = u 1 on3J u , (16b) 

an = t, ond&cr, (16c) 

where 38 is the domain of interest and the boundary subsets d£3 u and dSS a 
have essential and natural boundary conditions specified, respectively. These 
subsets satisfy dS§ u n = and dM^lJdM^ = dS§. 

2.3.2 Electric flow 

The electric flow problem is governed by the charge conservation equation with 
Dirichlet and Neumann boundary conditions: 

V • i = 0, m.38, (17a) 

^ = ^>, on d^, (17b) 

i ■ n = i, on dSBi, (17c) 

where d08$ n dSSi = and dW^ U dS§i = dSS. 

2.3.3 Mass Transport 

The mass-transport problem is governed by the continuity equation for vacan- 



cies, and appropriate initial and boundary conditions: 
3C 1 

— - = -V-j v --(C v -C?) X , m&, t>0, (18a) 
at t 

C V =C°, m.38, t = 0, (18b) 

C V =C V , on dm Cv , t > 0, (18c) 

j v -n=j, ond% v) t>0, (18d) 



where 83§c^ an d d&j are the boundary subsets on which concentration and 
flux boundary conditions are specified, respectively. Here, the boundary subsets 
satisfy d3$c v <~]d&j = and dSSc^ U — dSS. The effectiveness of vacancy 

sources and sinks is characterized by the relaxation time, r. The equilibrium 
vacancy concentration, C^ q , is defined by fJ, v \ci g — ^t? m GD- 
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Figure 2: Variation of e£ and across boundary regions. 



2.4 Grain-boundary migration 

The current location of a migrating grain boundary determines the value of 
the indicator function, x( x it) ( see Eq. ©); i.e. it determines whether vacancy 
sources/sinks are active and whether creep strain can accumulate at a given 
point, x 6 33. Furthermore, information about the location of the boundary 
is needed to calculate the value of the vacancy formation energy, e(, and the 
activation energy for diffusion, eij, everywhere in the domain of interest. These 
properties are assumed to vary linearly over the width of a boundary region as 
shown in Fig. |3 

The foregoing illustrates the influence of grain-boundary migration on mass 
transport. Due to the tight coupling between mass transport and mechanics, 
the migration of the grain boundary also affects the stress. In turn, mass transfer 
across the grain boundary causes one grain to grow at the expense of its neighbor 
and thus leads to grain-boundary migration. 

2.4.1 Thermodynamic driving forces 

Generally, interface migration in polycrystals is driven by the accompanying 
decrease in the free energy of the system. The thermodynamic driving force for 
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such a process, acting on a unit area of the interface, is thus defined as 



S3 



(19) 



p = -w> 



where S3 is the increase in the total Gibbs free energy of the system brought 
about by a motion of the interface, during which the interface sweeps through the 
volume SV. Neglecting triple junctions, the Gibbs free energy of a polycrystal 
can be expressed as the sum of two contributions: 



where g is the Gibbs free energy density as defined in and 7, is the energy 
per unit area of an interface, IV The term interface is used here to refer to free 
surfaces as well as grain boundaries, and the summation in (|2U[1 is over all such 
interfaces in the polycrystal. 

In situations where the Gibbs free energy density, g, suffers a decrease across a 
grain boundary, the total free energy of the system can be reduced if the grain 
with the smaller value of g (evaluated at the grain boundary) were to grow at the 
expense of its neighbor. Thus, a driving force acts on the grain boundary. For 
example, during recrystallization, annealed grains grow at the expense of cold- 
worked grains in which large dislocation densities lead to high values of g. The 
misorientation between two adjacent grains of an elastically anisotropic material 
subjected to a directional load causes one grain to store a smaller amount of 
strain energy per unit volume than its neighbor. This leads to strain- induced 
grain-boundary migration. Electromigration also causes atoms to jump across 
grain boundaries and thus leads to the migration of these boundaries [21] ■ 

Under the current formulation, the driving force discussed above can be char- 
acterized as follows: When a single atom hops across the grain boundary, it 
exchanges positions with a vacancy. Therefore, the volume change associated 
with this hopping event is SV = (1 — f)fl. The decrease in the Gibbs free en- 
ergy of the system is equal to the difference in the chemical potential of atoms 
across the boundary; i.e. — S3 = A^ a , or equivalently — S3 = — A/i„. Hence, 
from Eq. (15}, and assuming that the chemical potential gradient across the 
boundary is essentially linear [21], the driving force, p g , can be expressed as 



where n is the unit normal to the grain boundary and 2S g b is its width. 

From Eq. (|2(jp. it is clear that the total Gibbs free energy of a polycrystal can 
also be lowered by decreasing the total surface area of the grain boundaries in 
the system. When a curved interface moves away from its center of curvature, 




(20) 



Pg = 



(i-m 



-Aylt„ 



(- V/i„ ■ n)2S gb 

(i-m 



(21) 
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Figure 3: Driving force on a cylindrical grain boundary with a radius of curva- 
ture, R, due to surface tension, 7 5 b. 



sweeping through an increment of volume 5V, the free energy of the system 
increases by 

* s =^(i+£K (22) 

where 7$ is the (constant) specific surface energy of the interface and R\ and 
i? 2 are its principal radii of curvature. This expression can be traced back to 
Herring -26]. It follows that, in the current two-dimensional formulation, the 
driving force acting on a unit area of a (cylindrical) grain boundary, due to this 
effect, can be expressed as 

V, = , (23) 

where R is the radius of curvature and the negative sign indicates that p 1 
drives the boundary to migrate toward its center of curvature (see Fig- EJ - hi 
a polycrystal comprising only annealed grains, the action of p 1 leads to normal 
grain growth; i.e. the growth of large grains at the expense of smaller ones. By 
contrast, in the early stages of recrystallization, small annealed grains grow at 
the expense of the surrounding matrix of cold- worked material. In this case, the 
boundaries surrounding the annealed grains move away from their respective 
centers of curvature under the effect of p g . 



2.4.2 Kinetic law 

Assuming that grain-boundary migration takes place as a result of individual 
atoms hopping across the boundary, the migration velocity can be expressed in 
the following form: 

v n = Mp, (24) 
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where v n is in the direction of the local unit normal, n, which is assumed to point 
away from the boundary's center of curvature and M is the grain-boundary mo- 
bility. This classical result, derived by Turnbull [23] using absolute reaction rate 
theory, holds provided that pil <c kT, a condition which is always met in grain 
growth and recrystallization [27j. The results of molecular dynamics simula- 
tions of curvature-driven [2S| and strain-induced |2H] grain-boundary migration 
in bicrystals agree with Eq. 1)24(1. which can be modified as follows, to account 
for both types of driving forces discussed in Section Fi. 4. II 

V„ = MgPg + MyP^. (25) 



2.4.3 Grain-boundary mobility 

In Section 12.4.11 it was established that diffusion of atoms across a grain bound- 
ary due to the local gradient in the atomic chemical potential causes the bound- 
ary to migrate. The driving force for boundary migration in this case, denoted 
p g , is given by l|21l) . An expression for the corresponding mobility, M g , can be 
obtained by stipulating that, in this case, the migration velocity in the direction 
of the local unit normal, n, should be given by (see Porter and Easter ling 30 ) 

MgPg = - C?' a • n) fi - (j v ■ n) fn, (26) 

where j a and j v are, respectively, the local atomic and vacancy fluxes. Also 
recall that Q is the atomic volume and fQ is the volume of a vacancy. Since 
atoms and vacancies move by exchanging positions; i.e. j v = —j a , we have 

M gPg = (j v ■ n) (1 - f)n. (27) 

Finally, by combining (|27|l . (|21|l and l|15|) . we obtain 

M °- 2^fcT ■ (28) 

Since grain-boundary migration under the sole influence of p g involves transport 
across the boundary, and since, in the present treatment, mass transport is 
assumed to take place through the exchange of positions between atoms and 
vacancies, no lattice sites are transferred across the boundary in this process. 
By contrast, migration under the influence of p 7 consists of the transfer of lattice 
sites across the grain boundary |27|. via the process of atoms detaching from 
one grain and attaching to the one on the opposite side of the boundary. It is 
hence unreasonable to expect the mobilities, M g and M 7 , of these two distinctly 
different processes to be the same. 

As mentioned earlier, normal grain growth can be attributed to the action of 
the driving force p 7 . The fact that grain growth is a thermally activated pro- 
cess suggests an Arrhenius-type relationship between the mobility, M 7 , and the 
temperature: 

M 7 = M exp('^y (29) 



11 



where the pre-exponential factor, Mq, and the activation energy for grain- 
boundary migration, e m , are dependent on the misoricntation angle and the 
axis of rotation |3J [32] ■ 

Finally, by combining J2BJ), (|25|) . (|27[1 and l|29|l . we obtain the following expres- 
sion for the migration velocity: 

v n = (j v • n) (1 - f)n M cxp 'M. (30) 

3 Computational methods 

In this section, we focus on the level-set formulation of the grain-boundary mi- 
gration problem, and we present a detailed description of the computational 
techniques used in its implementation. To solve the coupled problem, the 
level-set formulation is integrated into the computational framework that was 
first introduced in Reference p^. This computational framework is based on 
an operator-split solution scheme and relies on the finite element method to 
solve the mechanics, mass-transport and electric-flow problems individually. 
While Reference J] addresses mainly the physics of stress-driven mass transport 
in polycrystals, it also contains details of the computational framework used 
therein. Here, these details are omitted; however, the operator-split solution 
scheme is outlined briefly in Section 13.21 to show how the level-set formulation 
is incorporated into the computational framework used. 



3.1 Level-set formulation 

Evolving interfaces can be tracked using the level set method, originally intro- 
duced by Osher and Sethian |16j . A comprehensive review of the method and the 
computational algorithms used in its implementation can be found in [331 134| . 

Consider an evolving grain boundary, T, which divides the domain of interest, 
SS 1 into two disjoint open subsets, 3$~ and S$ + . This situation is depicted in 
Fig. 2] The boundary can be parameterized with the aid of the scalar function 
(j)(x,t), defined on 88, provided that the following conditions are satisfied for 
all t > 0: 

cj)(x,t)<0, VseJ", (31a) 

cj>(x,t)=0, VxeT, (31b) 

<f>(x,t)>0, Vse (31c) 

The term level set refers to a set of points with a fixed value of (j>, i.e. an iso- 
contour of </>; the zero level set represents the grain boundary. Accordingly, the 
unit normal to a given level set can be defined locally as 

n + = TiSL (32) 
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Figure 4: Schematic of the domain of the grain-boundary migration problem. 



where || • || denotes the Euclidean norm. This expression can be evaluated at 
any point on the zero level set to obtain the local unit normal to the boundary, 
Tig". This definition implies that rig" always points into the SS + region. 

The evolution of the level-set field is governed by 

^ +J F„||V0||=O, (33) 

where F n (x,t) is the (scalar) local propagation velocity of the level set passing 
through point x. To track the motion of the grain boundary, we require that 

F n (x,t)=v n (x,t), Va;er, i > 0, (34) 

where v n is obtained at any point on the grain boundary from l(30|) . Although 
this requirement does not place any restrictions on the choice of F n away from 
r, the solution procedure is simplified greatly if F„ is an extensional velocity 
field; i.e. if 

VF n • n+ = 0. (35) 

This first-order partial differential equation can be solved for F n , at all points 
x G (3§\T), using 1)34(1 as a boundary condition. A simpler alternative strategy 
for constructing the extensional field is discussed in Section ?6 . 1 . II below. 

The level-set function, cp, can be initialized as the signed distance from T as 
follows: 

(j)(x, 0) = ^min \\x - p\\^j sign[(a; - p) ■ n+(p)]. (36) 

It is noted that this initial condition satisfies (|31|l . It also implies that initially, 

HV0H = 1, Vxe m, (37) 

Importantly, this desirable mathematical property of the level-set field is pre- 
served when the velocity field is extensional. This can be shown (see Jl]) by 
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noting that 



|||V^ = |(V0.V^) 

= 2V0 ■ ^ V0. (38) 
Combining 138|) and 1|33[) and assuming that </> and i 7 ^ are smooth, we obtain 
^||V0|| 2 - -2V</> ■ VF„||V0|| - ■ V||V0||F„. (39) 

Thus, if F n is extensional (i.e. V</>- VF n = 0) and </> is initially a signed-distance 
function (i.e. ||V0|| = 1 and V||V0|| = 0), we have 

^IIW||=0, (40) 
which implies that (|37|l holds for all t > 0. Hence, Eq. I|33|l reduces to 

f + ^ = o, («, 

which governs the evolution of <fi(x,t) from the initial condition 1)36(1 . 

Recall that the mechanics, mass-transport and electric-flow problems are solved 
using the finite element method, and that part of the coupling between these 
three sub-problems, on one hand, and the grain-boundary migration sub-problem, 
on the other, is through the indicator function, x( x i t)- I n this setting, the level- 
set update is computed at each node via a generalized trapezoidal rule: 

4>{xA,t n+ \) = (j){xA,t n ) ~ AtF n (x A ,t n+a ), (42) 

where xa is the position vector of node A and t n+a = at n+ \ + (1 — a)t n , with 
< a < 1. It is noted that when explicit time integration is used (a = 0), 
the update operation is trivial since v n , and hence F n , are determined from the 
known solution at t = t n and do not depend on (j)( x , t n +i)- It is also noted that 
no stabilization is required due to the use of the simplified level-set equation 
<|41|) in lieu of (|33|) . The advantages as well as the performance and stability 
characteristics of the resulting level-set formulation are studied in detail by 
Mourad et al. (17) . 

Additionally, \ is defined more precisely as follows: 

X {x,t) = H{8 gb -\(j>{x,t)\), (43) 
where H(-) is the Heaviside function. 
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Figure 5: Construction of the extensional velocity field by velocity projection. 



3.1.1 Velocity projection and field re-initialization 

The approach adopted here for constructing the extensional level-set propaga- 
tion velocity field, F n , is based on the notion that F n (x) — v n (p) for any x £ T, 
if p is such that 

1 1 as — p\\ = min \\x — p\\ . (44) 
per 

It is clear that the resulting velocity field satisfies (|35|l when 

( « - nHir (45 > 

as depicted in Fig. [S] This is not always the case however; for instance, if 
x G SS^ (see Fig. [SJ, then p is such that 

\\x-p\\= min ||x-p||. (46) 
pe(rna^) 

It follows that, for all x £ F n {x) = const., i.e. VF n = 0, which clearly 
satisfies l|35[) also. The above arguments apply in as well. 

Due to the accumulation of numerical error, the level-set field may develop 
perturbations; i.e. || V</>|| may deviate from unity in some regions within 3$. The 
field must be re-initialized to neutralize these perturbations and retain accuracy 
by maintaining ||V0|| = 1 (see Eqs. (jSHJEHJ)- Box ^ shows the algorithm 
used to implement the velocity projection scheme outlined above and the re- 
initialization scheme based on Eq. . 

It is noted that, to preserve the integrity of the solution, the re-initialization 
procedure should not change the current location of T. In other words, the re- 
initialized field should have the same zero level set as the original (perturbed) 
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field. The re-initialization scheme described here does not satisfy this condition 
strictly, i.e. it introduces error into the solution. Re-initialization should there- 
fore be used judiciously. This issue is examined in detail by Mourad et al. [*T7). 
An alternative re- initialization procedure designed to minimize this type of error 
is described in detail in [351 136] . It must be noted however that, in some cases 
including curvature-driven migration, quadratic convergence in L 2 is achieved 
with the level-set update formula 141fl . and importantly, this optimal conver- 
gence rate is preserved by the present re-initialization scheme (see ^| %Z7\ for 
details). 



FDR each node, A, DO 




SET D[A] = +00 




SET F[A] = 




ENDDD 




FDR each segment, Lq , of the zero level set 


DO 


FOR each node, A, (with position vector 


x A ) DO 


FIND point q 6 Lq such that 




\\xa -q\\ = min \\x A - q\\ 




q 6 La 




IF \\xa ~q\\ < \D[A]\ THEN 




COMPUTE v n (q) using Eq. (I3U1) 




COMPUTE n$(q) using Eq. (32) 




SET F[A] = v n (q) 




SET D[A] = \\x A - g|| sign[{x A - q) ■ n+(q)} 


END IF 




ENDDO 




ENDDO 




FOR each node^ A, DO 




IF x A £ {38+ U S%~ ) THEN 




RE-INITIALIZE = D[A] 




END IF 




ENDDO 





Box 1: Velocity projection and level-set field re-initialization algorithm. Here, 
the global arrays, €>, F and D hold the nodal values of 0, F n and \\x — p\\, 
respectively. A division of the zero level set which spans a single element is re- 
ferred to as a segment and is denoted by Lq. Also, since bilinear shape functions 
are used, an element can contain only one such segment. 

The algorithm in Box ^ was first introduced by Malladi et al. [3Hj and was pre- 
viously used by Garikipati and Rao It is adopted here for its simplicity, 
despite being relatively expensive — the complexity of the present algorithm is 
at best Q(riL x n np ), where is the number of elements intersected by the 
zero level set and n np is the total number of nodal points in the mesh (see [IT]). 
A more efficient algorithm, such as the fast marching method [331 BT)] with 
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0(n np log n np ) complexity, could be employed for velocity projection and level- 
set field re-initialization on Cartesian grids. 



3.1.2 Gradient smoothing 

The local migration velocity, v n , on the grain boundary is dependent on the 
local curvature, k = l/R, which is defined as follows: 

n = V-n+. (47) 

From ll.')2[l and (|37|) . it is clear that n + = V0, and the curvature can hence be 
expressed as 

where n sc i = 2 is the number of spatial dimensions. Since the value of <j> is 
updated at the finite element nodes, it is convenient to use the shape functions to 
evaluate the spatial derivatives in the above expression. However, since bilinear 
shape functions are used for simplicity and robustness, Eq. I|48[) cannot be used 
to evaluate n directly. To overcome this difficulty, we introduce a 'smoothed' 
normal vector field n + , weakly related to by 

J w ■ (n + - V</>) dV = 0, (49) 

m 

where w is an arbitrary weighting function. Equivalently, the nodal values 
of each component, fif ', of the smoothed normal vector can be obtained by 
minimizing the following discretized functional: 



E 



^ .N A nt(x A ) - ~^^ XA . 



E 

A=l 



dV, (50) 



where n e i is the number of elements in the model, SS e denotes an element 
domain, n en is the number of nodes per element and Na is the shape function 
associated with node A. This leads to a matrix equation of the form Md = f , 
where d is the global vector containing the nodal values of the component hf . 
The global mass matrix, M, and right-hand side vector, f, are obtained from 
the corresponding element arrays via the usual assembly process. The element 
arrays in this case are given by 



m AB = J N A N B dV, (51a) 
fl= [ ^E^^ s )dV. (51b) 



B=l 
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Finally, the curvature is evaluated as follows: 



7h en lL s d 

-EE 



dN A 



dxi 



n+{x A ). 



(52) 



A=l i=l 



Remark 2 Using this least-squares technique to smooth the stress field leads to 
a mixed problem in stress and displacement form (see Zienkiewicz and Taylor |41p 
By analogy, using this approach to compute the curvature is formally equivalent 
to a two-field mixed formulation for <fr and fi + . 

Remark 3 Similar techniques have been used previously to evaluate the cur- 
vature of an evolving interface; e.g. see Chessa and Belytschko |42|. 

3.2 Operator-split algorithm 

The coupled problem is solved using an operator-split algorithm. The sequence 
of operations carried out in one time step is shown in Box to illustrate how 
the level-set formulation is incorporated into this solution scheme. 



(1) CONSTRUCT the smooth level-set gradient field, n + 

(2) CONSTRUCT the extensional velocity field, F n , and 
RE-INITIALIZE the level-set field (see Box©. 

(3) PERFORM the level-set update using Eq. d42t . 

(4) SOLVE the electric-flow problem for the electric 
potential, ip . 

(5) REPEAT the following sequence: 

(a) SOLVE the mechanics problem for the 
displacements, u. 

(b) SOLVE the mass transport (composition) 
problem for the vacancy concentration, C v . 

UNTIL both the mechanics and mass transport 
problems have converged. 

(6) INCREMENT time and GOTO step (1). 



Box 2: The operator-split algorithm used to solve the coupled problem. 



4 Numerical examples 

Numerical results, obtained by solving the coupled initial and boundary value 
problem, are presented here with the aim of highlighting some of the advantages 
of the current approach. A 1 /im wide, 2.5 /im long segment of an aluminum 
interconnect line is modeled. The values of the material parameters used for Al 
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Table 1: Material properties of Al (thin film) used in the analysis 



Parameter 


Value 


Unit 


Elastic modulus, en 


184.7 


GPa 


Elastic modulus, C12 


95.15 


GPa 


Elastic modulus, C44 


44.7 


GPa 


Linear coefficient of thermal expansion, a 


24 x 1(T 6 


K- 1 


Electric resistivity, p 


4.2 x 10" 8 


Q,.m 


Apparent electric charge on a vacancy, q 


5.6077 x 1(T 19 


C 


Atomic volume in the absence of strain, Qo 


16.61 


A 3 


Vacancy-atom volume ratio, / 


0.8 




Vacancy formation energy in the bulk, e[ 


0.67 


eV 


Minimum value of el in boundary regions* 


0.5433 


eV 


Activation energy for diffusion in the bulk, e% 


1.47 


eV 


Minimum value of e% in boundary regions* 


1.1090 


eV 


Grain-boundary width, 2& gh 


0.198 


/im 


Diffusivity premultiplier^ , D VQ 


2.6 x 10 3 


m^s" 1 


Reduced grain-boundary mobility premultiplier* , Ao 


39.81 


m^s" 1 


Activation energy for grain-boundary migration, e m 


1.29 


eV 



* See Fig. 13 



The diffusivity is given by D v = D V(j exp(— ejJ/fcT). 
A = 7 9 i,M ; also see Eqs. |25J and ESJ- 



are given in Tabled The segment consists of two pure Al crystals separated by 
a £7 tilt grain boundary (38.2° misorientation about <111>). 

The line is assumed to operate at T — 373 K, with a reference temperature, 
T = 473 K. Rigid passivation material surrounding the line prevents vacancies 
from crossing the upper and lower boundaries of the domain (see Fig.EJi); i.e. 
the condition j v -n = holds at these boundaries. Periodic boundary conditions 
are imposed on the vacancy concentration at the left and right boundaries of 
the segment and an electrostatic potential difference, Aip = 0.0021 V, is applied 
between these two extremities. Vacancies drift along the electric field, E — 
— V^, pointing to the right. 

In the first example, the migrating grain boundary consists, initially, of two 
straight (planar) sections, which are joined by a circular (cylindrical) section. 
The straight sections form 45° angles with the upper and lower boundaries 
of the domain as shown in Fig. Since the vacancy formation energy, el, 
is low inside grain-boundary regions where vacancy sources are also present, 
vacancies accumulate in such regions resulting in high values of the local vacancy 
concentration, C v . The location of the grain boundary is revealed locally by the 
maximum- valued contour of C v . 

Figure |S] shows the evolution of the vacancy concentration contours in the line. 
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Vacancy Concentration 
[cm -3 ] 

■ -8.34E-01 

- 6.80E+00 

1.44E+01 

_ 2.21E+01 

2.97E+01 

3.73E+01 

_ 4.50E+01 

* 5.26E+01 



Time = 30G.0 sec 




Figure 6: The evolution of the vacancy concentration contours in the intercon- 
nect line due to the motion of the grain boundary (first example), (a) t — 2.0 sec. 
(b) t = 300.0 sec. (c) t = 900.0 sec. 
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Flux Magnitude 
[cm _2 .s _1 ] 




O.OOE+OO 
1.48E-01 
2.97E-01 
4.46E-01 
5.94E-01 



7.43E-01 



8.91E-01 



1.04E+00 



Time = 300.0 sec 



Figure 7: Contours of the magnitude of the vacancy flux, \\j v \\, in the line (first 
example) . 



It is clear that, during the first 300 seconds, the central cylindrical section of 
the boundary migrates to the right; i.e. toward its center of curvature, while 
the planar sections are relatively less mobile (Fig. EJd) . It must be emphasized 
that the interaction between mechanics, mass transport, electric effects and 
grain-boundary motion is accounted for, and the driving force for boundary 
migration, due to stress-driven diffusion and electromigration, is included in the 
calculations. However, its effect is overshadowed by the dominant driving force 
due to the curvature of the boundary. 

After 900 seconds, the curvature is approximately the same everywhere on the 
grain boundary (Fig.Efc). At this stage, a steady state prevails and the boundary 
continues to travel toward the right without undergoing any further changes in 
shape. It is noted that the level-set calculations remain stable. The boundary 
remains smooth and does not develop any spurious cusps or ripples. The 45° 
equilibrium angles between the grain boundary and the sidewalls are maintained 
as the solution progresses. 

A contour plot of the magnitude of the vacancy flux at t = 300 sec is shown 
in Fig. It is clear that a strong vacancy flux exists in the vicinity of the 
grain boundary. A vector plot of the vacancy-flux field in the neighborhood 
of the grain boundary is shown in Fig. [S] and although the field is complicated 
in this neighborhood, it can be seen that vacancies tend to drift in the same 
direction as the migrating boundary, thus preventing the appearance of a 'trail' 
of vacancies or other oscillations in the numerical solution in the boundary's 
wake. It is important to note that the flux field evolves continually as the 
boundary migrates. It is also notable that the formulation captured this aspect 
of the coupling between mass transport and grain-boundary migration without 
additional terms being added to the expression of the vacancy flux to account, 
specifically, for the effect of moving boundaries. 

The second example is concerned with the case where a curved grain boundary 
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Figure 8: Vector plot of the vacancy flux, j v , in the neighborhood of the grain 
boundary at t = 300 sec (first example). Also see Fig. [7| for magnitude of j v . 



evolves into a planar configuration to reduce the free energy of the system. Here, 
the equilibrium angles at the grain boundary-sidewall intersections are set to 
90° and the initial geometry of the boundary is different and less regular than 
in the first example, ft is clear from Fig. [51 which shows the evolution of the 
vacancy concentration contours in this case, that the boundary flattens as the 
solution progresses. It is noted that the equilibrium angles are also preserved 
in this case. It is also noted that in this case, perturbations in the solution 
lead to the formation of spurious 'shadow' regions , @t~ ; see Fig. |SJ and 
re-initialization becomes necessary in these regions to maintain accuracy and to 
preserve the contact angles. 

Although the numerical stability characteristics of the level-set formulation — 
and those of the operator-split algorithm — are not examined in a formal setting, 
no spatial or temporal oscillations are observed in the numerical solution of the 
example problems presented, as long as the time-step size is within the CFL 
limit; i.e. At < h/F n , where h is the mesh parameter. A detailed study of the 
performance and stability characteristics of the present level-set formulation is 
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Vacancy Concentration 
[cm -3 ] 

_ -2.37E-02 

W- 2.71E-01 

1 5.66E-01 

H 8.61E-01 

^ 1.16E+00 

1.45E+00 

_ 1.75E+00 

* 2.04E+00 



Time = 2.00 sec 





Time = 2.00 sec 



Figure 9: The evolution of the vacancy concentration contours in the in- 
terconnect line due to the motion of the grain boundary (second example), 
(a) t = 2.0 sec. (b) t = 200.0 sec. (c) t = 1500.0 sec. 
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presented by Mourad et al. ^7j. Details regarding the advantages, applications 
and numerical stability characteristics of operator-split schemes can be found 
in |431 1441 14"5] and references therein. 

5 Summary and Concluding Remarks 

In this paper, we present the following contributions. 

• A computational formulation capable of capturing the full coupling between 
grain-boundary motion and other microscale phenomena that take place in 
pure polycrystalline materials. The coupled continuum formulation pre- 
sented here was developed to model stress-driven self-diffusion and elec- 
tromigration in polycrystals while accounting fully for the interaction be- 
tween these mass transfer processes and the motion of grain boundaries. 
The formulation accounts for two distinct thermodynamic driving forces 
acting on a grain boundary; one due to the boundary's own curvature and 
another engendered by mass transfer across the boundary via stress-driven 
self-diffusion and electromigration. 

• A simplified lev el- set formulation which does not require spatial stabilization. 
The level set method has previously been used to pose grain-boundary 
migration as a time-dependent field problem governed by a pure advec- 
tion equation. Standard numerical schemes resort to spatial stabilization 
techniques (e.g. upwinding schemes, Galerkin/Least-Squares) to attenu- 
ate the spurious oscillations known to appear in the numerical solution 
of equations of this type. Here, on the other hand, the level-set equation 
is reduced, using the mathematical properties of signed distance func- 
tions and extensional velocity fields, to a simpler form which obviates the 
need for these stabilization techniques. This leads to a remarkably sim- 
ple explicit scheme for advancing the solution in time. The algorithm 
used to construct the extensional velocity field and to simultaneously re- 
initialize the level-set field — is presented. We also provide the description 
of an L 2 -projection technique used to compute the curvature of the grain 
boundary. 

The numerical examples presented indicate that the strong coupling in the prob- 
lem is captured adequately and that the numerical implementation allows the 
solution of the coupled initial and boundary value problem to be advanced in 
time in a stable fashion to obtain physically meaningful results. 
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